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We propose a phase space method to propagate a quantum wavepacket driven by a strong external field. 
The method employs the so-called biorthogonal von Neumann basis recently introduced for the calcula- 
tion of the energy eigenstates of time-independent quantum systems [A. Shimshovitz and D.J. Tannor, 
arXiv:1201. 2299vl]. While the individual elements in this basis set are time- independent, a small subset 
is chosen in a time-dependent manner to adapt to the evolution of the wavepacket in phase space. We demon- 
strate the accuracy and efficiency of the present propagation method by calculating the electronic wavepacket 
in a one-dimensional soft-core atom interacting with a superposition of an intense, few-cycle, near-infrared 
laser pulse and an attosecond extreme-ultraviolet laser pulse. 



With the emergence of attosecond laser technology^ 
there is the fascinating prospect of observing and con- 
trolling the correlated dynamics of multiple electrons on 
its natural time scale of ten to one hundred attoseconds. 2 
In order to unravel the complex and sometimes counter- 
intuitive^ quantum dynamics from the experimental 
data, and to develop theories that reproduce the essence 
of the dynamics^, an accurate and efficient numerical 
method to simulate the multi-electron wavepacket dy- 
namics is indispensable. 

However, accurate simulation of the electronic dynam- 
ics in a high-intensity laser field is a challenging task: the 
electronic wavepacket is dispersed by the laser field over a 
wide region of coordinate space while retaining high mo- 
mentum near the atomic nuclei. Straightforward repre- 
sentation of the wavefunction on a equally-spaced coordi- 
nate grid [i.e., a Fourier grid (FG)] requires a large range 
with a small interval between points. Simulation on such 
a large grid quickly becomes prohibitive as the number 
of degrees of freedom (DOF) increases. Even with so- 
phisticated treatments such as multi-configuration time- 
dependent Hartree-Fock (MCTDHF),— simulation of ion- 
ization dynamics have been limited to small systems such 
as the helium atom and the hydrogen molecule^— 

In this article, we present a new approach to solv- 
ing the time-dependent Schrodinger equation (TDSE), 
based on a phase space perspective. The resultant prop- 
agation method is simple, accurate, stable, and effi- 
cient. As a first demonstration, we simulate the elec- 
tronic wavepacket of a one-dimensional (ID) model atom 
in the combined fields of a high-intensity near-infrared 
(NIR) laser pulse and an attosecond extreme-ultraviolet 
(XUV) laser pulse. 

In our approach, we utilize the localized nature of 
phase space Gaussians to prune the basis. However, a 
basis set of phase-space localized states is perforce non- 
orthogonali^, and this has created difficulties in previ- 
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ous attempts^— to truncate the basis to cover only 
the phase space region of dynamics. This long stand- 
ing problem was solved recently, by the so-called peri- 
odic von Neumann (pvN) basisi 20 ' 21 The pvN basis is 
generated from a set of Gaussians whose centers form 
a finite lattice in phase space. By imposing periodic 
boundary conditions on these Gaussians, the pvN ba- 
sis becomes formally equivalent to the simple and accu- 
rate FG representation^ Then, its biorthogonal basis, 
the biorthogonal von Neumann (bvN) basis, can be used 
for a compact representation of a quantum state. This 
framework was successfully applied to the calculation of 
quantum energy eigenstates^ Here, the framework is 
extended for solving the TDSE. Our strategy is to keep 
the individual elements in the basis set time-independent 
but to truncate the bvN basis in a time-dependent man- 
ner. This avoids the problem of a moving basis^ that 
can become over-complete as time elapses and can be- 
come unstable. In contrast, here we obtain a stable set 
of linear ordinary differential equations (ODEs) for the 
expansion coefficients in the truncated bvN basis. 

Before explaining the pvN and bvN bases, we re- 
view the formalism of the FG basis to establish nota- 
tion. We write the Fourier pseudo-spectral and spectral 
bases^-^as {|0 m )} m =i,...,jv and {\(j) m )} m =i,...,N , respec- 
tively. These bases are orthonormal and span the same 
A^-D Hilbert space denoted here as W, resolving the iden- 
tity in T-L as 

N N 

In = J2 \<MiM = E l<U<<U (1) 
m—1 m—1 

The bases {|0 m )} an d {|0m)} are localized at the 
FG points {a; m }"m=i,...,iV in the position space and 
{Pm}m=i,....N in the momentum space, respectively. For 
any quantum state \^) G %, {6 m \^) = {x m \fy)\J~Kx and 
(0 TO |\I r ) = (p m \^)y/Ap, where Ax and Ap are the grid 
intervals in position and momentum spaces, respectively. 
The pvN basis^&^i {\9j}}j=i,...,N is defined as 

N 

\~9j) = E \O m ){x m \gj)V^, (2) 

m—1 
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where {\gj)}j = i n are the phase space Gaussians, 



\9j) 



2h mj 



(3) 



whose centers {((jj,^)}^! at constitute a finite lattice 

in the phase space with the unit cell of area 2irh and the 
momentum-to-position aspect ratio fvy. Note that the 
number of Gaussians N is the same as the number of FG 
points used. 

The bvN basis^i {\bj)}j = i at is defined to be 

biorthogonal (dual) to the pvN basis, i.e., (bi\9j) = 
This gives 



N 
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where S 1 is the inverse of the overlap matrix = 



(gi\gj) = AxJ2 m =i(9i\ x ™)( x m\gj) of thepvN basis. The 
pvN and bvN bases span the same Hilbert space H as 
{|^m)} an d {|0m)}, resolving the identity as 



N N 

In = Eiffel = £&>& 

3=1 3=1 



(5) 



Thus, any wavepacket € % can be represented as 



N 



i*w) = Ei6i>(&-i*w>- 



(6) 



Due to the localized nature of the Gaussians {|<?j)}, the 
magnitude of (gj\^f(t)} can be extremely small if the 
corresponding classical system can not reach the phase 
space region around (qj,pj). Defining a set A such that 
|(<?j-|^(i)}| is negligible if j ^ A, we can approximate the 
wavepacket by a subset of the bvN basis as 



E fe>< 



(7) 



where Cj(t) := (gj\^(t)}. Note that the set A of active in- 
dices can be changed in time in order to keep the number 
Na of elements in A small at all time. 

By substituting cq. Q to the TDSE, we obtain a set 
of linear ODEs for the active bvN coefficients {cj}j^_A, 



dc 3 

dt 



= -i£ E (^hibiwrn^it), (8) 



where SI 1 is the inverse of the overlap matrix flji = 
(bj\bi) = (S~ 1 )ji of the bvN basis, and H(t) is the Hamil- 
tonian operator of the system. The overlap and Hamilto- 
nian matrix elements in eq. ((SJ can be computed simply 
via the representations in {|# m )} and {|^m)}' 21 ' 23 The 
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FIG. 1. Vector potentials, -Anir(£) and Ax_vv(t), of the NIR 
and XUV laser pulses applied to the model ID atom. 



matrix SI -1 is Hermitian positive-definite, and the el- 
ements {(bi\H (t)\b m )} constitute an Hermitian matrix. 
Therefore, the product of these matrices yields all real 
eigenvalues^, and eq. ((SJ) can be solved stably by many 
standard numerical algorithms. 

To demonstrate the accuracy and efficiency of the 
present method, we solve eq. © for the electronic 
wavepacket of a ID atom in the combined field of NIR 
and XUV laser pulses. The Hamiltonian of this system 
is given as 



H(t)=H + V(t), 
where Hq is the field-free Hamiltonian expressed as 



(9) 



Qe 2 



4tt€ \/x 2 + a 2 



(10) 



Here fi — 1 a.u. is the electron mass, e = — 1 a.u. is 
the electron charge, — Qe = 1 a.u. is the charge of the 
atomic nucleus, a = 1 a.u. is the soft-core parameter, 
and e = 1/Att a.u. is the electric constant. The laser- 
electron coupling V{t) is, in the velocity gauge, 



V(t) 



;[4m(t) + A xvv (t)] P , 



(11) 



where Amir^) and Axw{t) are the vector potentials of 
the NIR and XUV laser pulses, respectively. We used an 
NIR pulse of wavelength 800 nm, peak intensity 5 x 10 13 
W/cm 2 , and duration 1.5 cycles (4 fs, FWHM of intensity 
profile). The XUV pulse had wavelength 15 nm, peak 
intensity 1 x 10 12 W/cm 2 , and duration 5.0 cycles (250 
as). The peak of the XUV pulse was delayed from that 
of the NIR pulse by 0.25 NIR cycles, as shown in Fig. [T] 
The initial state was chosen as the ground state 
of Hq with energy eigenvalue —0.66978 a.u., and the 
wavepacket was propagated from the turn-on of the NIR 
laser pulse at t m j n to its end at t max by the short-iterative 
Arnoldi algorith m 27 ' 28 in a 6D Krylov space with a con- 
stant time-step of At = 0.0379 a.u. We divided the time 
span from t m - m to i max into 8 time segments and changed 
the active set A from one segment to the next. The 
number of FG points was JV = 4096, and they were dis- 
tributed over —750 a.u. < x < 750 a.u. and —8.58 a.u. 
< p < 8.58 a.u. This phase space rectangle was divided 
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FIG. 2. Snapshots of {|cj| 2 }j6^ shown by the ellipses located 
at the Gaussian centers {(lj,Pj)}j€A- The colors of the el- 
lipses indicate the magnitude of \cj\ 2 according to the scale 
above the figure. The sequence of dark blue dots represent 
the simple- man trajectories for direct ionization; the light blue 
dots represent the rescattered simple-man trajectories. The 
dark blue + marks represent the simple-man trajectories ab- 
sorbing one XUV photon in the presence of the NIR field. 
The snapshots were taken at (a) t = —2.06, (b) t = 0.69, and 
(c) t = 2.06 in units of NIR cycles. These times are indicated 
by the green x marks in Fig. [1] 
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FIG. 3. Comparison of the photoelectron momentum distri- 
butions obtained with the reduced bvN basis (blue solid line) 
and full bvN basis (red dashed line). The momentum distri- 
bution from a simulation without the XUV pulse (using the 
full bvN basis) is also shown (gray solid line). The vertical 
dashed lines indicate the cut-offs of the direct (Nl and NT') 
and rescattered (N2 and N2') photoelectrons, as well as the 
NIR-streaked single-XUV-photon ionization peaks (XI and 
XI'), estimated by the simple-man model. 



into 64 x 64 cells of aspect ratio fry = 1.14 x 10 a.u., 
each of which contained a Gaussian center (qj,pj). 

In Fig. [2[ snapshots of |cj(i)| 2 are shown by the color 
scale of the ellipses located at the active Gaussian cen- 
ters {(qj,Pj)}jeA- The outer rectangular boundary of 
each panel indicates the phase space area correspond- 
ing to the Hilbert space H spanned respectively by the 
Fourier spectral and pseudospectral bases as well as the 
full pvN and the full bvN bases. The wavepacket, initially 
concentrated at the atomic core [Fig. [Ha)], is ionized by 
the NIR and XUV laser pulses and spreads into parts of 
the first (x > and p > 0) and third (x < and p < 0) 
quadrants [Fig. [2fb,c)], but a large area is never accessed. 
This can be intuitively expected from the classical me- 
chanics, and indeed we see that the wavepacket closely 
follows the so-called simple-man trajectories^ (dots and 
+ marks in Fig. [2]) which obey the classical Hamilto- 
nian of the same form as eq. © with Hq replaced by 
p 2 /2/j,. In fact, we chose A so that the active phase space 
domain contains these simple-man trajectories (with an 
additional margin) . 

In Fig. [21 we compare the photoelectron momentum 
distributions obtained using the reduced basis of Fig. [2] 
and the full bvN basis. The excellent agreement between 
the two results indicates that the present method not 
only preserves the qualitative features — cut-offs of the 
direct and rescattered NIR photoelectrons, and the NIR- 
streaked single-XUV-photon ionization peaks — but also 
has quantitative accuracy. 

The accuracy and the size of the active bvN basis 
set are expected to be inversely related. However, it 
is not straightforward to determine the minimal size of 
the active set that maintains a given accuracy. To seek 
an upper bound to such an optimal basis size, we car- 
ried out simulations with different active sets generated 
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FIG. 4. The error e as a function of (N A )/N (black x marks). 
The horizontal error bars indicate the range of N A (t)/N in 
tmin < t < t m ax- The data marked by the red filled circle is 
from the simulation shown in Figs. [2] and [3] 



by changing the margins around the simple-man tra- 
jectories. The error of a simulation was measured by 
e := |(* rcduc (t,„a X )|* full (ima X ))-l|, where |* rcduc (W)> 
and I* (tmax)) are the final states calculated using a 
reduced and the full bvN bases, respectively. Figure 
H] shows the dependence of this error on the basis set 
size. As the size Nj, is time-dependent, and as the 
computational cost of wavepackct propagation depends 
on JV4 quadratically, we characterized N A by its root 

mean-square, (NX, 



max L mmji 

as well as its minimum and maximum values. It can 
be seen that the bvN basis can be compressed down to 
(Na)/N = 0.14 or possibly further while maintaining 
the accuracy level at e = 4 x 1CP 10 , and at least to 
(N A ) /N = 0.08 for e = 6 x 10~ 8 . 

The computational cost of the present method is dom- 
inated by the multiplication of the matrix Gj m := 
(-i/fyY^ieAi^^jl&ilHfylbm) and the vector c m 
(j, m £ A) which scales as 0(N A ). There is also initial 
overhead that scales as 0(N 3 ) originating from the com- 
putation of S m j and (S~ 1 )ji (m,j = 1, ...,N; I £ A). 

Our method can be extended straightforwardly to sys- 
tems with multiple DOF.— Defining the number d of 
DOF, the total number M of FG points increases ex- 
ponentially as M = 0(N d ). In contrast, the number M A 
of active bvN basis states for atomic and molecular sys- 
tems in a laser field is expected to scale more slowly than 
M since the phase space coordinates are generally corre- 
lated for such systemSi 4 ' 29 i 30 The cost per time step of the 
present method goes as 0(M A ). The initial overhead in 
calculating S m j and (S~ 1 )ji stays at 0(N 3 ) because the 
multi-D Gaussians can be factored into ID Gaussians^. 
The dominant part in the overhead is now the calculation 
of the matrix elements for the potential energy operator 
in the reduced bvN basis, and this goes as 0(MM A ). 
For sufficiently large d, the present method has the po- 
tential to perform better than the popular alternate- 
direction Crank-Nicolson schem o 31 ' 32 [which costs O(M) 
per step] and the FG split-operator mctho d 24 i 25 [which 
costs 0(M log 2 N) per step]. The bvN basis may also be 
used for the one-body orbitals in MCTDHF as an alterna- 



tive to the static scaling of the position coordinate d 8 ' 33 ^ 4 
used previously^ 

In summary, we presented a new method to solve the 
TDSE based on the bvN basis. Although the basis is 
time- independent, the active subset is chosen in a time- 
dependent manner. As a first demonstration, we calcu- 
lated the electronic wavepacket of a ID atom in the com- 
bined fields of intense NIR and attosecond XUV laser 
pulses. This example demonstrates the high accuracy 
and efficiency of the method. We are currently working 
to extend the method to 3D with the aim of ultimately 
applying it to multi-electron systems in intense and ul- 
trashort laser pulses. 
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